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1.  INTRODUCTION 


The  Naval  Surface  Warfare  Center  (NSWC),  in  collaboration  with  Germany,  has 
developed  a  computer  simulation  program,  (DYSMAS),  to  model  both  the  linear  and 
nonlinear  behavior  of  ship  structures  under  shock  loading.  This  program  is  targeted  towards 
computer  simulation  of  U.S  Navy  ship  shock  trials.  The  failure  process  must  be  modeled  in 
an  accurate  and  reliable  way  in  order  to  meet  the  objectives  of  this  program.  One  effort  to 
model  this  nonlinear  failure  process  is  to  implement  damage  constitutive  equations,  like 
Gurson’s  void  model,  into  the  program. 

Shell  elements  are  used  for  modeling  since  the  major  structure  in  a  ship  is  a  shell. 
The  formulation  of  this  shell  element  must  include  damage  constitutive  equations.  The 
element  must  accurately  reflect  the  elastic-plastic  transition  in  a  dynamic  environment. 
Therefore,  the  overall  objective  of  this  research  is  to  develop  a  shell  element  which  can 
include  the  damage  constitutive  equations. 

The  first  phase  of  the  research  developed  a  shell  element  which  can  include 
Gurson’s  void  model  at  the  next  phase.  Since  Gurson’s  void  model  uses  hydrostatic  pressure 
(mean  stress)  and  deviatoric  stress,  these  values  must  be  available.  The  second  phase  is  to 
implement  the  void  model  into  the  shell  formulation  developed  in  the  first  phase.  The  final 
phase  is  to  incorporate  the  entire  module  into  the  DYSMAS  program.  Each  phase  requires 
extensive  verification  of  the  developed  shell  element. 

The  first  phase  was  completed  in  September,  1997.  All  applicable  portions  of  the 
report  covering  the  first  phase  are  repeated  in  this  report,  since  some  of  the  procedures 
outlined  there  had  to  be  modified  to  implement  Gurson’s  void  model  [1].  This  report  covers 
the  work  of  both  the  first  and  second  phases. 
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2.  FINITE  ELEMENT  FORMULATION 


2.1  Geometry 

A  point  in  a  shell  structure  can  be  expressed  by  a  vector  sum  of  two  vectors.  The 
first  vector  is  a  position  vector  from  the  origin  of  the  coordinate  system  to  a  point  on  a 
reference  surface  of  the  shell  element.  The  second  vector  is  a  position  vector  from  this 
reference  surface  to  the  point  under  consideration.  The  surface  that  spans  the  center  of  the 
transverse  axis  is  used  as  the  reference  surface  in  this  formulation,  although  any  surface  would 
suffice.  The  first  vector  terminate;  i  the  reference  surface  directly  below  the  point  in 
question.  The  second  vector  is  then  the  normal  from  the  reference  surface  that  intersects  the 
desired  point.  Figure  1  shows  this  relationship. 


Figure  1.  Element  Cross-section. 
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Two  shape  functions  are  used  to  describe  a  position  in  the  element;  A*  is  the  two 
dimensional  shape  function  in  the  £-r|  plane,  and  Hk  is  the  one  dimensional  shape  function 
along  the  £  axis,  where  (£,r|,£)  describes  a  point  in  the  natural  coordinate  system.  A  generic 
point  of  a  shell  may  now  be  described  in  terms  of  the  position  vectors  of  the  nodes  and  the 
shape  functions: 

*,(5,1.0  =  E  Nk(^,r\)xjk  +  E  *‘(5,1)# ‘(0  Vi  (>  =  1,2,3)  (1) 

k= 1  k=\ 

k  k 

where  x.  is  the  position  vector  of  node  k  in  the  reference  surface;  V3i  is  the  unit  vector  at  the 
node  k. ;  and  n  is  the  number  of  nodes  per  element.  In  the  present  formulation,  a  four-node 
shell  element  is  considered.  The  unit  vector  V3i  is  defined  as 


„  (*  ,r  -  (*, ) 


bottom 


v;,  = 


(x,y°p  -  <x‘y 


k\bottom 


(2) 


where  top  and  bottom  indicate  the  top  and  bottom  surfaces  of  the  shell,  and  ||  ||  denotes  the 
Euclidean  norm.  The  one-dimensional  shape  function  Hk  is  expressed  as 


H\ 0  = 


1(1  +  CXI  -0(1+0 

4  4 


k\bottom 


(3) 


in  which  £  indicates  the  location  of  the  reference  surface  and  varied  from  -1  to  1  (  £  =  0 
denotes  the  mid-surface).  The  two-dimensional  shape  function  A*  is  expressed  as 

JV  =  7(1 -5X1-1) 

N2  -  7(1  +5)(1  “l) 

=  7(1  +5)(i  +i)  <4) 

4 

A4  =  1(1-  0(1  +T1) 
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2.2  Displacement 

The  displacement  field  in  a  shell  can  be  written  as 


=  (•  -  i,2,3)<5) 

£=1  £=1 


in  which  w;  is  the  displacement  along  the  x,  axis,  uf  is  the  nodal  displacement  at  the  node  k, 

k  k  kkk 

and  unit  vectors  Vu  and  V2i  lie  along  the  reference  surface.  Vu,V2i  and  V3i  are  perpendicular 
to  one  another.  0*,  0^  and©*,  are  rotational  degrees  of  freedom  along  the  unit  vectors 

k  k  k 

Vy,V2i  and  V3i,  respectively.  The  right-hand  rule  is  assumed  for  the  positive  direction  of 
each  rotation.  Figure  2  illustrates  the  relationship  among  these  vectors. 


Figure  2.  Displacement  Vector  Orientation. 
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0^  and  0j;  are  the  bending  rotations,  while  0^  is  the  drilling  rotational  degree  of 
freedom.  The  role  of  0^  can  be  clarified  by  considering  a  flat  plate  parallel  to  the  XjX2  plane. 
Now  Equation  (5)  can  be  rewritten  as 

ui  =  u™d  +x30;  (i  =  1,2,3)  (6) 

where  and  u2  are  the  inplane  displacements,  and  u3  is  the  transverse  displacement.  The 
superscript  mid  indicates  the  mid-plane  of  the  plate.  Equation  (6)  demonstrates  that  the 
transverse  displacement  is  not  constant  through  the  plate  thickness  (i.e.  along  the  x3  axis). 
Therefore,  the  transverse  normal  strain  is  included  in  this  formulation,  along  with  the 
transverse  shear  strains. 

In  the  implementation,  the  unit  direction  vectors  are  fixed  to  the  first  node  listed 
in  the  element,  where  the  nodes  are  ordered  in  a  counter-clockwise  direction,  as  shown  in 
Figure  3.  V3  is  normal  to  the  plane  formed  by  nodes  1,  2,  and  4.  V2  lies  along  the  line 
connecting  nodes  1  and  4.  V1  is  normal  to  the  plane  formed  by  V2  and  V3 .  This  means  that 
Vj  will  not,  in  general,  lie  along  the  line  connecting  nodes  1  and  2.  For  each  node,  the 
following  information  is  stored  at  each  time  step:  displacement  {d^  dy,  dP  0x,  0^,  0Z),  velocity 
(v„  vy  6x,  dy  0z),  and  acceleration  ( a ^  dy,  a.,  Qx,  6y,  0Z).  Using  a  four-node  element,  each 
node  has  six  degrees  of  freedom  (dof),  and  each  element  has  24  dof  s  in  displacement, 
velocity,  and  acceleration. 


Figure  3.  Node  Number  and  Unit  Direction 
Vector  Scheme. 
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2.3  Coordinate  Transformation 

Combining  the  three  unit  direction  vectors  into  matrix  Tp  provides  the  rotation 
transformation  matrix,  or  matrix  of  direction  cosines,  as  shown  in  Equation  (7).  T'1  is  used 
to  transform  the  nodal  angular  displacements  from  the  global  coordinate  system  to  local 
coordinates,  as  shown  in  Equation  (8),  where  k  is  the  node  number.  Once  the  internal  force 
vector  is  generated  in  local,  Tp  is  used  to  transform  the  moments  back  into  global  coordinates. 
This  procedure  is  discussed  later  in  section  2.7. 


V  V  V 

v\  ’ 1  k3 


1(3*3) 


(7) 


10  0  0 
0  10  0 
0  0  10 
0  0  0 

0  0  0 
0  0  0 


(k  =  1 ,2,3,4) 


(8) 


The  strain  transformation  matrix,  T,  is  used  to  transform  calculated  strain  from  the 
global  coordinate  system  to  local  coordinates.  Transforming  the  resulting  stress  from  local 
coordinates  to  the  global  coordinate  system  would  normally  require  using  T  but  since  T  is 
orthogonal,  T'1  =  TT,  where  TT  is  the  transpose  of  T.  This  property  negates  the  requirement 
to  invert  a  six  by  six  matrix.  For  a  detailed  derivation  of  these  transformations,  refer  to  Cook 
[2],  The  strain  transformation  matrix  is  explicitly  defined  in  Equation  (9),  where  Vy  is  the 
cosine  of  the  direction  vector  V;  in  the  xy  direction.  Both  of  these  transformation  matrices 
are  calculated  once  per  element  per  time  step. 
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2.4  Strain  Displacement  Relation 

The  six  components  of  the  strain  tensor  are  computed  from  Equation  (5)  by  taking 
its  derivative  with  respect  to  the  x,  axis.  In  matrix  form,  the  result  for  a  four-node  element 
is: 

{e}  =  [S]{4  (10) 

where 

{e}  =  jcjj  e22  e33  y12  y23  y13j  (11) 


[B]  =  [Bl  B2  B*  BA]  (12) 

{< d }  =  {d1  d2  d 3  d4}T  (13) 
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The  detailed  expression  for  [Bk  ]  is 


M  - 


dNk 

dxl 

0  0 

kjrk 

Si  V2l 

kjrk 

Si  Vn 

kjrk 

Si  V31 

0 

o 

dx2 

kjrk 

Si  ^22 

kjrk 

Si  Vn 

kjrk 

Si  V32 

0 

0  *!i 

-SiV^ 

sX 13 

kjrk 

S\  V33 

dNk 

Sx2 

0 

ax, 

-sXi-sXi 

Si^+gXi 

kjrk  kjrk 

Si  V31  +Si  V32 

0 

dNk  dNk 
axj  axj 

kjrk  kjrk 

Sk  *22  ~Sk  *23 

kjrk  kjrk 

Sk  ^12 +Sk  ^13 

kjrk  kjrk 

Sk  V31  +Sk  V33 

dNk 

3x3 

0 

ax. 

-sXi  -sX3 

kjrk  kjrk 

S3  Vn  +Si  Vn 

kjrk  kjrk 

S3  ^31  +Sl  ^33 

(14) 


in  which 


gk  =  dNk_jjk  +tfkdH 

*v_ 


3X; 


dx; 


(15) 


The  vector  {cf  }  is  defined  as 


{d-}  -  {,*  «,*  ej  e‘  e^} 


(16) 


where  w;  is  the  displacement  along  the  x,  direction  at  node  k,  and  0,  is  the  rotational 
displacement  about  the  x,  axis  at  node  k.  The  matrix  [Bk  ]  must  be  calculated  for  each 
integration  point.  The  stress  tensor,  discussed  later,  is  also  stored  as  a  vector. 
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2.5  Jacobian  Matrix 

Computing  the  derivatives 


—  and  —  requires  the  Jacobian  matrix,  defined  as 

dxt.  dxi 


M- 


X2,C 

xu 

x„ 

x. 

l.T] 

2,r\ 

3,r| 

.X« 

X2A 

X3,C 

where 


%-±2£x,*±2£H'rj  (,  =  1.2.3) 

35  h  35  '  t,  35  ' 


JL*'E  ^f~x,  *  E  ~Hkvi  (,=  1.2.3) 
3r]  jk=i  5r|  it=i  dr^ 


dx,. 


9-  =  T;xk^v£  0  =  14,3) 

d(  t=i  dt. 


(17) 

(18) 

(19) 

(20) 


[i?]  is  defined  as  the  inverse  of  the  Jacobian  ( J Then  the  required  partial 
derivatives  are  defined  as 


dNk  D  dNk  D  BNk 

=  /?., - +R „ -  (/  =  1,2,3) 


dx. 


as 


,2 


at] 


dHt-R.,™l  (,  =  1.2.3) 


dx; 


,5 


as 


(21) 

(22) 
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Both  the  Jacobian  and  its  inverse  are  calculated  at  each  integration  point  for  each 
time  step.  The  inverse  of  the  three  by  three  matrix  J  is  explicitly  calculated  (no  loops). 


2.6  Stress-Strain  Relationship 

The  strain  calculated  in  Equation  (10)  is  in  global  coordinates,  and  is  transformed 
to  the  local  coordinate  system  using 


{6Zoca/}  lT]{€g/oia/}  (23) 

where  [7]  is  defined  in  Equation  (9).  Stress  is  calculated  from  the  strain  (local  coordinates) 
using  the  plane-strain  formulas: 


■v = 

E  K+eJ 


a  i  = 

y 


l-v: 


a  „/  =  Ee, 


V/  =  Gy 'xy 


Xy'z<  =  KGyyz 


Xxu  =  K  Gyx 


(24) 


where  K  is  the  shear  correction  factor,  E  is  the  elastic  modulus,  G  is  the  shear  modulus,  and 
v  is  Poisson’s  ratio.  The  resulting  stresses  are  in  the  local  coordinate  system  and  are 
converted  to  the  global  coordinate  system  with 

{»}  =  MV)  (25) 


\o\  =  lo  o  o  y  y  y  \ 

l  >  \  x  y  z  1  xy  1  yz  *  xz) 


(26) 


where  T  is  from  Equation  (9). 
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2.7  Internal  Force,  Mass,  and  Assembly  into  System  Matrix 

The  stress  resulting  from  Equation  (25)  is  then  converted  to  an  internal  force  vector 
and  summed  over  all  integration  points  using 


1  i  i 


nx  ny  nz 


\f internal)  =  f  f  f  [BY M <%<**<%  =  £  E  £  [J3f{o}w.W  wj/|  (27> 

ixix  -i  ,=l  ;=l  fc=1 


where  |  J  |  is  the  determinant  of  the  Jacobian  matrix,  nx,  ny,  and  nz  are  the  number  of 
integration  points  in  the  xh  x2,  and  x3  directions,  respectively,  and  wh  wp  and  wk,  are  the  Gauss 
weights  for  those  directions.  The  resulting  force  vector,  must  have  its  moments 
transformed  back  to  the  global  coordinate  system  using 


1  0  0  0  0 

0  10  0  0 

0  0  10  0 

0  0  0 

0  0  0  [rj 


0  0  0 


( k 


1  >2,3,4) 


(28) 


A  lumped  mass  method  is  used  for  the  element  mass  matrix.  The  matrix  is 
diagonal,  with  each  diagonal  element  the  same: 


iw  0  0 

e 

0  m  0 

e 

0  0  m 

e 


o 

o 

0 


0  0  0...  7W 

e 


(29) 


11 


where 


nx  ny  nz 


EEE 


„  w:» 

e 


(30) 


w 


with  w  =  4  in  this  four-node  element.  Using  a  diagonal  mass  matrix  greatly  simplifies  time 
integration,  since  inverting  it  is  trivial  and  takes  little  computation  time.  The  internal  force 
vector  and  element  mass  vector  (each  24  by  1)  are  then  assembled  into  the  corresponding 
system  vectors,  each  of  which  is  a  column  vector  having  one  row  for  each  dof. 
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3.  EXPLICIT  TIME  INTEGRATION 


The  use  of  internal  force  vectors  and  explicit  time  integration  negates  the  need  to 
explicitly  form  the  element  and  system  stifihess  matrices.  The  acceleration  vector  is 
computed  from 

{0}'  =  [M]-'({F4'-{F„}')  (31) 

where  {f/}  is  the  system  acceleration  vector,  [M]  is  the  system  mass  matrix,  is  the 
system  external  force  vector,  and  superscript  t  denotes  the  time  step.  Velocity  and 
displacement  are  then  quickly  found  using 

M'**  =  {£/}'^ *{(/}'  ,32) 

{C/}f+Ar  =  {C/}f  +  {f/}  2  At 

The  boundary  conditions  are  applied,  and  the  data  is  written  to  an  output  file,  along  with  the 
stresses,  strains,  yield  stresses,  and  void  contents  for  each  integration  point. 
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4.  DAMAGE  CONSTITUTIVE  EQUATIONS 


4. 1  Gurson’s  Void  Model 

Yielding  and  plastic  deformation  in  the  element  follows  the  model  proposed  by 
Gurson  [3]  for  a  symmetric  deformations  around  a  spherical  void.  The  yielding  condition  is 


0  = 

(  \ 

JL 

,  % 

2 

+  2<7j/cosh 

'_3q# 

9  n 

^  0  / 

-(l*?/2)=  0 


(33) 


where  /  is  the  current  porosity,  p  is  the  hydrostatic  stress,  q  is  the  effective  stress,  and  a0  is 
the  current  tensile  yield  stress.  The  constants  qh  q2,  and  q3  are  constants  introduced  by 
Tvergaard  [4]  in  order  to  provide  a  better  match  with  numerical  studies.  Aravas  [5]  provides 
a  detailed  explanation  of  implementing  this  model  in  a  static  finite  element  algorithm.  The 
procedure  used  here  is  essentially  the  same,  with  minor  modifications  due  to  the  different 
element  formulation.  The  stress  is  then  transformed  to  local  coordinates,  and  the  internal 
force  vector  is  computed  as  described  above  in  section  2.6.  One  major  benefit  of  using  this 
model  in  Equation  (33)  is  that  if /is  initially  0,  the  yielding  criteria  surface  is  identical  to  the 
von  Mises  yield  condition.  Addessio,  et  al.,  [6]  discussed  the  importance  of  void  growth  and 
its  effects  on  plastic  deformation  in  high  strain  rate  conditions,  especially  on  the  pressure 
wave  propagation  through  the  material.  Therefore,  this  model  is  well  suited  to  the  application 
of  ship-shock  test  simulation. 

After  calculating  the  strain  tensor  using  Equations  (24),  any  previous  plastic  strain 
is  subtracted  using 


{ee}  =  {etotal}  -  {eP} 


(34) 


The  stress  is  then  calculated  using  the  components  of  {ee}  in  Equation  (24).  Using  these 
values,  the  hydrostatic  stress  and  effective  stress  are  calculated  as  follows: 
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(35) 


p  =  —  (a  +a  +o\ 

^  3  l  »  yy  &  I 

{5}  =  {o}+/>{5}  {5}  =  {l  1  1  0  0  0}r  (36) 

q  =  +2(s*  +s*  +s*))  (37) 


At  this  point,  O  is  calculated  from  Equation  (33).  If  0  is  greater  than  0,  indicating  plastic 
deformation,  iteration  is  required  to  determine  the  new  porosity  and  change  in  plastic  strain. 
The  predictor-corrector  method  used  by  Addessio,  et.  al.  [6],  is  also  used  here.  Using  the 
values  of/?  and  q  previously  calculated  as  a  first  guess,  correction  factors  are  calculated  using 

{Cf}  =  Wl{Al}  (38) 


where 


{Al}  = 


(39) 


K3*_ 

dp 


—  +KAe 

dg 
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c2® 

dp 2 


-3G— 

dg 


— -3GAen 
dp  p 


e2® 

dg2 


(40) 


(41) 
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and 
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where  fN  is  the  volume  fraction  of  void  nucleating  particles  and  eN  and  Sn  are  the  mean  and 
standard  deviation  of  a  normal  distribution  of  nucleation  strain,  as  suggested  by  Chu  and 
Needleman  [7],  and  utilized  by  Aravas  [5].  These  values  are  included  in  the  input  file  element 
property  matrix.  Then  the  effective  plastic  strain,  and  void  content  are  updated  with 


AeP 


-PAep+gAeq 

(i-/K 


(P 

ti-At 


ef  +Aep 


(47) 


where  ef  is  the  effective  plastic  strain  from  the  previous  time  step.  At  this  point,  the  change 
in  yield  stress  due  to  strain  hardening  is  calculated,  which  is  discussed  in  section  4.2.  If  either 
a  or  P  is  greater  than  a  predetermined  tolerance,  the  process  iterates  beginning  with  Equation 
(39).  In  practice,  the  above  procedure  is  rarely  required  more  than  once  due  to  the  critical 
time  step  limitation  of  explicit  time  integration.  The  tolerance  used  is  1.0  x  10"4 .  If  the 
critical  time  step  size  is  exceeded,  the  resulting  stress  will  diverge  towards  infinity,  and  the 
iterative  procedure  outlined  above  will  not  converge.  If  the  number  of  iterations  exceeds  a 
predetermined  value,  the  element  algorithm  forces  the  program  to  terminate  with  an 
appropriate  error  message.  Using  a  calculation  time  step  size  less  than  the  critical  time  step 
size  will  allow  convergence  of  the  plastic  deformation  algorithm.  During  the  iteration 
process,  the  updated  values  of  effective  plastic  strain,  Ae  A  e  f  and  {e^}  are  retained  and 

r  7 

used  as  starting  points  for  the  next  iteration. 

4.2  Strain  Hardening 

Modeling  the  nonlinear  elasto-plastic  behavior  of  the  material  used  is  simplified  by 
constructing  a  piece-wise  linear  version  of  the  stress-strain  plot.  As  written,  the  element  will 
support  up  to  15  separate  line  segments  to  model  the  elasto-plastic  region.  Input 
requirements  are  further  simplified  by  allowing  the  input  of  the  tangent  modulus,  ET,  directly 
from  the  tensile-test  stress-strain  plot,  and  the  upper  strain  limit  of  the  region  being  modeled. 
The  new  yield  stress  is  calculated  with 
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(48) 


°o+A'  =  °o  +  E  Er,(e>  -  ef-i)  +  Erjft  ~  ej- 1) 

j=i 


where  A>  is  the  yield  stress  for  this  time  step,  a0is  the  original  yield  stress,  ef  is  the 
effective  plastic  strain  for  the  time  step  under  consideration,  and  e(is  the  upper  strain  limit  of 
the  z'th  linear  segment.  e0is  the  original  yield  strain,  and  is  calculated  by  the  program  as 
simply  — .  Figure  4  illustrates  an  example  calculation. 

E 


Figure  4.  Calculating  a  New  Yield  Stress. 

The  effective  plastic  strain  from  Equation  (47)  is  frmnd  to  he  within  the  second 
plastic  segment.  The  new  yield  stress  is 

°0  =  a0  +ETj{ei  ~  6o)  +  ET2^\^t* Kt  ~  ei)  (49) 

The  algorithm  calculates  the  new  yield  stress  as  a  function  of  the  cumulative  effective  plastic 
strain  and  the  original  yield  stress  at  each  iteration. 
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Figure  4  also  shows  how  the  linear  model  for  the  elastic-plastic  region  is 
constructed.  For  most  materials,  three  or  four  segments  are  sufficient  to  capture  the  strain¬ 
hardening  behavior.  The  algorithm  allows  the  user  to  input  a  negative  tangent  modulus, 
which  is  required  to  accurately  model  steels,  necking,  and  failure  (discussed  in  section  4.3). 
The  material’s  behavior  in  compression  is  assumed  to  be  identical  to  its  behavior  in  tension. 
However,  with  only  a  minor  modification  to  the  program,  a  separate  piece-wise  linear  model 
for  compression  could  be  included.  The  maximum  number  of  linear  segments  has  been 
arbitrarily  set  to  fifteen,  but  this  can  also  be  readily  increased  to  accommodate  a  more 
complex  model. 

A  material  that  is  ideally  plastic  after  yielding  is  modeled  by  entering  one  linear 
segment  of  slope  zero  and  an  upper-limit  strain  higher  than  the  expected  strain.  If  the 
effective  plastic  strain  of  the  element  exceeds  the  upper-limit  strain  of  the  last  linear  segment, 
the  yield  stress  will  remain  constant  at  the  value  of  the  endpoint  of  the  last  segment.  No 
special  indications  are  provided  when  this  happens,  so  the  stress-strain  model  should  be 
defined  for  effective  plastic  strains  well  beyond  those  expected  in  the  structure.  An  important 
property  that  is  built  into  the  strain  hardening  and  plastic  deformation  algorithm  is  the 
retention  of  any  previous  plastic  strain.  If  the  external  force  is  removed,  the  plastic  strain  will 
remain.  As  load  is  removed,  the  stress-strain  plot  will  decrease  at  slope  E  from  the  new  yield 
stress,  and  the  new  yield  stress  will  remain  in  effect.  If  load  is  reapplied,  plastic  yielding  will 
not  occur  again  until  the  effective  stress  exceeds  the  new  yield  stress.  For  impact  modeling, 
this  property  is  vital  to  properly  predict  behavior  as  shock  waves  propagate  back  and  forth 
through  the  material. 

4.3  Failure 

The  strain  hardening  algorithm  used  in  this  element  formulation  also  allows  simple 
modeling  of  element  Mure.  Since  the  strain  hardening  algorithm  will  accept  negative  values 
for  tangent  modulus,  the  stress-strain  plot  can  be  continued  down  to  a  stress  of  nearly  zero 
in  any  manner  that  is  appropriate  to  reflect  the  material’s  failure  behavior.  A  final  yield  stress 
of  zero  is  not  currently  supported,  and  will  cause  a  divide  by  zero  error,  but  can  easily  be 
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allowed  if  desired  later.  (The  error  will  only  occur  if  the  element  is  strained  to  that  point;  the 
element  will  function  normally  at  any  strain  below  the  point  where  the  yield  stress  becomes 
zero).  Since  the  material  does  not  simply  disappear  when  failure  occurs,  a  small  amount  of 
yield  stress  should  be  retained  to  represent  the  momentum  and  stiffness  of  the  remaining 
structure.  Further  study  is  required  to  determine  an  appropriate  lower  limit  for  this 
representative  yield  stress. 

4.4  Storage  Concerns 

The  damage  constitutive  equation  model  requires  several  key  data  from  the 
previous  time  step  in  order  to  complete  the  calculations.  The  array  used  to  pass  the  stress  and 
strain  values  back  to  the  time  integration  routine  is  also  used  to  store  the  previous  values  of 
yield  stress,  porosity,  effective  plastic  strain,  and  the  plastic  strain  tensor  (stored  as  a  vector). 
The  elastic  stress  and  strain  is  simply  written  to  an  output  file,  and  is  not  required  for  future 
calculations.  A  column  vector  of  length  21  is  required  for  each  integration  point:  six  elastic 
stress  components,  six  elastic  strain  components,  six  plastic  strain  components,  current  yield 
stress,  current  porosity,  and  current  effective  plastic  strain.  The  program  structure  currently 
supports  up  to  a  total  of  eight  integration  points  per  element.  The  basic  formulation  is 
designed  for  one  in  the  x,  and  x2  directions:  the  “center”  of  the  element’s  plane  surface.  The 
number  of  integration  points  in  the  x3  direction  is  user-defined  in  the  input  file.  Currently,  all 
elements  must  have  the  same  number  of  integration  points. 

Several  constants  must  be  provided  in  the  input  file  for  this  model:  the  initial  void 
content,  the  constants  ql,  q2,  and  q3,  the  piecewise  linear  stress-strain  relationship,  and  the 
distribution  constants  for  void  nucleation.  The  current  implementation  does  not  utilize  default 
values,  which  are  therefore  provided  by  the  preprocessor. 
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5.  HOURGLASS  MODE  CONTROL 


5.1  Hourglass  Effects 

In  the  standard  formulation,  the  element  is  integrated  fully  along  the  x3  direction, 
but  is  under-integrated  in  the  xrx2  plane.  It  is  necessary  to  prevent  shear-locking  of  the 
elements,  which  results  in  a  structure  that  is  too  stiff.  If  any  component  of  the  applied  force 
lies  parallel  to  the  x,  -  x2  plane  of  any  element  in  the  structure,  the  zero  energy  mode  of  that 
element  is  excited  enough  to  adversely  affect  the  results.  This  excitation  is  quickly  propagate 
to  the  surrounding  elements,  and  eventually  to  the  entire  structure.  These  modes  are  referred 
to  as  “Hourglass  Modes”  due  to  the  shape  the  structure  takes  on  after  complete  propagation 
of  this  distortion.  The  impact  of  exciting  these  modes  is  to  completely  distort  the  structure 
in  an  unnatural  manner,  which  renders  any  output  -data  useless.  In  cases  where  the  applied 
force  does  not  sufficiently  excite  a  zero  energy  mode,  the  analysis  gives  reliable  results.  When 
modeling  a  general  curved  structure,  however,  it  is  almost  impossible  to  avoid  exciting  this 
mode.  Figure  5  clearly  illustrates  the  phenomenon  using  the  pinched  cylinder  verification 
problem  from  section  6.3.  In  this  case,  the  hourglass  mode  control  built  into  the  algorithm 
was  disabled,  and  distortion  in  the  structure  is  clearly  visible. 


Figure  5.  Pinched  Cylinder  without  Hourglass  Mode  Control. 
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5.2  Method  of  Control 

Belytschko,  et  al.,  [8]  proposed  an  efficient  means  of  controlling  the  hourglass 
modes  of  a  similar  element.  The  method  described  uses  a  portion  of  a  stiffness  matrix 
generated  by  full  integration  in  all  directions  to  modify  the  stiffness  matrix  generated  by 
under-integration.  Although  this  formulation  does  not  use  a  stiffness  matrix,  a  similar 
approach  is  just  as  effective  in  controlling  these  modes. 

In  this  formulation,  only  one  hourglass  mode  presents  a  problem,  as  shown  in 
Figure  5.  Rather  than  fully  integrating  in  all  three  directions,  the  element  is  fully  integrated 
in  the  x,  -  x2  plane,  but  under-integrated  in  the  x3  direction.  The  relative  location  of  the 
integration  points  used  for  hourglass  control  is  shown  in  Figure  7.  The  procedure  described 
in  section  2  is  followed  to  generate  an  internal  for  .e  vector  related  to  these  four  integration 
points. 


Figure  6.  Hourglass  Mode  Control  Integration  Points. 

The  algorithm  for  calculating  strain,  stress,  and  then  force  for  the  hourglass  mode  control 
integration  points  is  identical  to  the  algorithm  used  to  produce  the  main  internal  force  vector, 
with  the  exception  of  plastic  strain.  No  plastic  strain  is  subtracted  from  the  strain  resulting 
from  Equation  (23),  and  the  damage  constitutive  equations  described  in  the  previous  section 
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are  not  utilized.  The  internal  forces  generated  from  the  two  integration  schemes  are  treated 
like  the  stillness  matrices  in  Belytschko,  et  al.  [7],  The  new  force  vector  is 


W -{/£“}♦  *{r.4  <S0) 

where 

{M  =  <S1> 


and 

ft ^ 

h  -  T  <52) 

The  variable  h  is  used  here  in  Equations  (50)  and  (51)  instead  of  the  e  used  in 
Belytschko,  et  al.  [7]  to  avoid  confusion  with  the  multiple  strains  required  in  this  formulation. 
The  variables  used  to  calculate  h  are  the  element  thickness  (t)  and  the  element’s  surface  area 
(A).  The  effect  of  r  follows  that  described  in  Belytschko,  et  al.  [7],  and  is  set  to  0.05.  The 
range  of  values  for  r  that  effectively  controls  the  hourglass  modes,  but  does  not  greatly  effect 
the  overall  element  stiffness  is  roughly  0.046  to  0.057  (determined  experimentally).  Since  the 
elements  are  in  arbitrary  orientation  in  3-D  space,  the  area  calculation  is  computed  as  the  sum 
of  the  area  of  the  two  triangles  formed  by  dividing  the  element  at  the  diagonal  between  nodes 
1  and  3: 


A 


(53) 


where 


Di  =  (*i  -*2)(x3  +  CVi  -T2)03  ~y2)  +  (zi  "z2>(z3  -z2> 

d2  =  (*!  -*4)(x3  -*4) +  (Ti  -y4)(y3  -y4)+(zi  ~z4)(z3  ~z4) 


(54) 
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(55) 


h  =  \I{X4-X3j  +{y4~yif  +{Z4~Z3j 

h  =  \J[xi  -  x4f +  {y\  ~y^f +  [zi  ~  z4f 

and  (x*  ,  yk ,  zk )  is  the  location  of  node  k  in  the  global  coordinate  system.  For  these 
calculations,  the  element  is  assumed  to  be  flat  (no  curvature  along  either  the  Xj  or  x2 
directions). 

5.3  Effectiveness  and  Impact  of  Hourglass  Mode  Control 

The  hourglass  mode  control  described  was  applied  to  the  same  pinched  cylinder 
problem  shown  in  Figure  5.  The  results  are  shown  in  Figure  7,  which  clearly  demonstrates 
the  effectiveness  of  the  this  method  in  controlling  the  zero  energy  modes.  Both  Figure  5  and 
Figure  7  use  the  same  structure,  applied  force,  calculation  time  step,  and  time  step  displayed. 


Figure  7.  Pinched  Cylinder  with  Hourglass  Mode  Control. 
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The  additional  error  introduced  by  including  hourglass  control  is  sufficiently  small 
that  it  does  not  alter  the  results  significantly.  All  verification  problems  analyzed  in  the 
following  section  were  completed  with  hourglass  control  enabled.  The  current 
implementation  uses  hourglass  control  consistently,  but  allows  easy  modification  to  make 
hourglass  control  a  user-defined  option. 

One  possible  modification  is  to  apply  the  hourglass  control  force  every  two  or  three 
time  steps,  vice  every  step.  Another  possibility  is  to  apply  hourglass  control  to  selected 
elements,  vice  every  element.  Both  modifications  would  reduce  processing  time,  but  require 
further  study  to  determine  if  these  changes  would  provide  an  adequate  reduction  in  the 
hourglass  modes. 


6.  VERIFICATION  EXAMPLES 
6. 1  Transformation  Matrix  Verification 

A  cantilever  plate  is  subjected  to  a  tip  load  of  0.4  N  distributed  evenly  along  its 
width.  The  elastic  modulus  is  200  GPa,  the  density  is  7850  kg/m3,  Poisson’s  ratio  is  0.29,  and 
the  yield  stress  is  50  kPa.  The  plate  is  2m  x  2m  x  10  cm  thick.  The  problem  is  solved  in  four 
orientations:  parallel  to  the  x-y  plane,  parallel  to  the  y-z  plane,  parallel  to  the  x-z  plane,  and 
skew  to  all  three  planes.  Two  integration  points  through  the  thickness  are  used.  Beam¬ 
bending  theory  gives  a  bending  stress  at  the  top  integration  point  of  10.3923  kPa.  All  results 
are  summarized  in  Table  I,  following  the  descriptions  of  each  run. 

Figure  8  shows  the  orientation  for  the  first  run.  Nodes  1,2,  and  3  are  clamped,  and 
the  force  is  applied  in  the  -z  direction  at  nodes  7,8,  and  9.  Figures  9  and  10  show  the 
displacement  of  the  center  tip  node  and  the  stress  at  the  top  integration  point  of  element  #1, 
respectively. 


Figure  8.  Cantilever  Plate  in  the  x-y  Plane. 
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x10'5  Node  #8 


Figure  9.  First  Run  Displacement  History. 
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Figure  10.  First  Run  Bending  Stress  History. 
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Figure  1 1  shows  the  orientation  for  the  second  run.  Nodes  1,2,  and  3  are  clamped, 
and  the  force  is  applied  in  the  x  direction  at  nodes  7,8,  and  9.  Figures  12  and  13  show  the 
displacement  of  the  center  tip  node  and  the  stress  at  the  top  integration  point  of  element  #1, 
respectively. 


Figure  11.  Cantelever  Plate  in  the  y-z  Plane. 
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Figure  12.  Second  Run  Displacement.  Figure  13.  Second  Run  Bending  Stress. 
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Figure  14  shows  the  orientation  for  the  third  run.  Nodes  1,2,  and  3  are  clamped, 
and  the  force  is  applied  in  the  -y  direction  at  nodes  7,8,  and  9.  Figures  15  and  16  show  the 
displacement  of  the  center  tip  node  and  the  stress  at  the  top  integration  point  of  element  #1, 
respectively. 


Figure  14.  Cantilever  Plate  in  the  x-z  Plane. 
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Figure  15.  Third  Run  Displacement. 
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Figure  16.  Third  Run  Bending  Stress. 
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Figure  17  shows  the  orientation  for  the  third  run.  Nodes  1,2 ,  and  3  are  clamped, 
and  the  force  is  applied  in  the  -y  direction  at  nodes  7,8,  and  9.  Figures  18  and  19  show  the 
displacement  of  the  center  tip  node  and  the  stress  at  the  top  integration  point  of  element  #1, 
respectively. 


Figure  17.  Cantilever  Plate  in  Skew  Orientation. 
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Table  I  summarizes  the  results  for  the  four  runs  described  above.  The  displacement 
and  stress  values  are  the  means,  calculated  by  dividing  the  peak  value  by  two.  The  difference 
between  run  number  four  and  the  other  runs  is  due  to  the  accumulation  round-off  error  in 
each  component  of  the  applied  force  and  location  of  the  nodes.  The  applied  force  had  a 
magnitude  of  0.3996  N,  vice  the  0.4  N  used  in  the  other  problems.  Note  that  the  analytic 
solution  assumes  no  stress  across  the  width  of  the  plate.  The  force  was  applied  as  a  linear 
ramp  from  0.0  N  to  0.4  N  over  the  time  interval  0.0  seconds  to  0.1  seconds,  then  held 
constant  at  0.4  N,  to  prevent  shock  effects  from  influencing  the  results.  This  could  not, 
however,  remove  all  dynamic  affects,  and  is  the  main  source  of  error  between  the  analytic  and 
finite  element  solutions.  The  magnitude  of  the  applied  force  was  chosen  to  keep  the  effective 
stress  throughout  the  structure  well  below  the  yield  stress  of  the  material,  restricting  the 
response  to  the  elastic  region. 


Run  Number 

Orientation 

TABLE  I 

Tip  Displacement  (m) 

Bending  Stress  (kPa) 

1 

x-y  Plane 

27.6840x10"* 

9.839 

2 

y-z  Plane 

27.6840x10-* 

9.839 

3 

x-z  Plane 

27.6840x1 O'6 

9.839 

4 

Skew 

27.6300x10^ 

9.608 

Analytic 

Static 

29.3x10-* 

10.392 

The  results  above  demonstrate  that  the  transformation  matrices  used  in  this 
formulation  are  effective  in  converting  the  displacements,  strains,  and  stresses  between  the 
global  and  local  coordinate  systems. 
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6.2  Elastic  Plate 

A  plate  clamped  on  all  four  sides  is  subjected  to  a  concentrated  force  at  its  center. 
The  elastic  modulus  is  10  msi,  the  density  is  0. 1647  slugs/in3,  and  Poisson’s  ratio  is  0.2.  The 
dimensions  of  the  plate  are  10  in  x  10  in  x  0. 1  in  thick.  The  yield  stress  is  set  high  enough  to 
ensure  a  completely  elastic  response.  Two  integration  points  through  the  thickness  are  used. 
The  applied  force  is  40  lbf.  The  finite  element  mesh  uses  symmetry  to  model  one  quarter  of 
the  plate.  Both  a  2x2  (4  element)  and  4x4  (16  element)  mesh  are  used  in  the  finite  element 
analysis.  Figure  20  shows  the  structure  for  the  4  element  mesh. 


Figure  20.  Clamped  Plate  4  Element  Mesh. 


Nodes  1,2,3 ,4,  and  7  are  clamped.  Nodes  6,  8,  and  9  have  symmetry  conditions 
applied.  The  force  is  applied  as  a  step  of  magnitude  -10  lbf  at  node  9.  The  results  for  both 
meshes  and  the  analytic  solution  are  shown  in  Table  DL 


TABL 

Analysis  Type 

EH 

Center  Node  Peak 

Displacement  (in) 

2x2  FE  Mesh  (Dynamic) 

-4.74xl0'2 

4x4  FE  Mesh  (Dynamic) 

-4.94xl0‘2 

Analytic  (Static  times  2) 

-4. 90x1  O'2 
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6.3  Pinched  Cylinder 

An  open-ended  cylinder  of  radius  5.0  in.,  length  10.35  in.,  and  thickness  0.094  in. 
is  subjected  to  a  pinching  load  of  100  lbf.  The  elastic  modulus  is  10.5  msi,  Poisson’s  ratio  is 
0.3125,  and  the  density  is  3.125xl0‘3  slugs/in2.  The  load  is  applied  as  a  step  function 
beginning  at  time  t  =  0  seconds.  Figure  21  illustrates  the  problem. 


100  lbf 


5.175 

5.175 

—  *  • 
i 

f 

100  lbf 


(All  Dimensions  in  Inches) 
Figure  21.  Pinched  Cylinder  Problem. 


Using  symmetry,  the  problem  was  reduced  to  a  one-eighth  section  of  the  cylinder. 
The  16  element  mesh  used  is  shown  in  Figure  22.  The  deflection  in  the  Z-direction  of  node 
number  25  is  shown  in  Figure  23,  which  corresponds  to  the  radial  contraction  of  the  cylinder. 
The  dynamic  value  should  be  twice  the  analytic  static  value.  Inextensional  shell  theory  gives 
a  static  radial  contraction  of  0. 1 1 17  in.  The  maximum  radial  contraction  of  the  model  is 
0. 1995  in.,  which  translates  to  a  static  radial  contraction  of  0.09975  in.  Using  a  256  element 
mesh  (16  by  16),  the  maximum  radial  contraction  was  0.2207  in.,  for  a  static  contraction  of 
0.1104  in. 
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Figure  22.  Pinched  Cylinder  Radial  Contraction. 
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6.4  Elastic-Plastic  Plate  without  Void  Nucleation 

A  structural  steel  plate  clamped  on  all  four  sides  is  subjected  to  a  shock  pressure 
in  the  form  of  a  step  with  amplitude  60  kPa.  The  elastic  modulus  is  200  GPa,  the  density  is 
7850  kg/m3,  Poisson’s  ratio  is  0.29,  and  the  yield  stress  is  250  MPa.  The  stress-strain 
relationship  will  model  the  results  of  a  typical  tensile  test.  The  plate  dimensions  are  9  m  x  9 
m  x  1  cm  thick.  A  9  element  mesh  (3  x  3)  will  be  used  to  model  one  quarter  of  the  plate.  The 
finite  element  mesh  is  shown  in  Figure  24.  Nodes  1, 2,  3, 4,  5,  9,  and  13  are  clamped.  The 
remaining  boundary  nodes  have  appropriate  symmetry  conditions  applied. 


Figure  24.  Finite  Element  Mesh  for  Elastic-Plastic  Plate. 

The  stress-strain  curve  for  the  lower  strain  region  is  shown  in  Figure  25.  A  larger 
region  of  the  stress-strain  curve  is  shown  in  Figure  26.  This  problem  is  designed  so  as  not 
to  exceed  the  maximum  stress  of  the  strain-hardening  region. 


35 


x  i q8  Lower  Region  of  Stress-Strain  Curve 


Figure  26.  Structural  Steel  Stress-Strain  Curve. 


Figure  25.  Structural  Steel  Stress-Strain  Relationship. 
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0  0.01  0.02  0.03  0.04  0.05 

Effective  Strain 

Figure  27.  Stress-Strain  at  Top 
Integration  Point. 


0  0.1  0.2  0.3  0.4  0.5 

Time  (sec) 


Figure  28.  Stress  History  at  Top 
Integration  Point. 


0  0.1  0.2  0.3  0.4  0.5 

Time  (sec) 

Figure  29.  Displacement  of  Center  Node. 


Figures  27  through  29  show  the  results  of  the  finite  element  analysis.  The 
effectiveness  of  the  strain  hardening  implemenation  is  shown  in  Figures  27  and  28.  With 
microvoid  effects  disabled,  the  stress-strain  curve  follows  the  tensile  test  curve  exactly.  Note 
in  Figure  27  that  the  stress  follows  the  elastic  modulus  down  as  stress  is  relieved  after  the 
peak  displacement. 
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6.5  Elastic-Plastic  Plate  with  Void  Growth  and  Nucleation 

The  problem  of  section  6.4  is  repeated  with  the  microvoid  effects  included.  The 
constants  for  Equation  (33)  are:  q1  =  1.5,  q2  =  1.0,  and  q3  =  2.25  (qx2).  The  initial  porosity 
is  0.0.  The  density  of  nucleating  particles  (fN)  is  4%,  the  mean  strain  for  nucleation  (eN)  is  0.3 
with  a  standard  deviation  (Sn)  of  0. 1.  These  are  the  values  recommended  by  Tvergaard  [4] 
and  Aravas  [5],  All  other  parameters  are  identical  to  the  problem  described  in  section  6.4. 


x  io9  Element  #3  Integration  Point  #1 
3  r - t - ; - i— 


2.5 


0  0.05  0.1  0.15  0.2 


Effective  Strain 

Figure  30.  Stress-Strain  Plot  at  Top 
Integration  Point. 


x  io*  Element  #3  Integration  Point  #4 


Fig;  31.  Stress- Strain  Plot  at  Bottom 
Intc  .aon  Point. 


Figure  33.  Center  Node  Displacement.  Figure  32.  Void  Nucleation  and  Growth. 
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Figures  31  and  32  show  the  effect  of  microvoids  on  effective  stress.  Figure  31 
shows  the  stress-strain  relationship  in  the  top  integration  point  of  element  #3.  This  fiber  of 
the  element  is  in  tension,  and  voids  quickly  nucleate  and  grow  during  plastic  deformation,  as 
shown  in  Figure  34.  Figuer  32  shows  the  stress-strain  relationship  at  the  bottom  integration 
point  of  the  same  node.  Since  this  fiber  is  in  compression,  there  is  no  void  nucleation,  and 
the  peak  stress  is  significantly  less  than  in  the  fiber  in  tension.  Comparing  Figure  33,  center 
node  displacement,  with  the  center  node  displacement  where  void  nucleation  is  disabled 
(Figure  30),  illustrates  that  the  peak  displacement  response  is  virtually  unaffected  by  the 
introduction  of  microvoids.  However,  as  Figure  32  shows,  the  element  experienced 
significant  additional  plastic  deformation  and  strain  hardening.  The  top  fiber  of  this  element 
is  much  closer  to  failure  than  the  corresponding  fiber  in  the  problem  where  microvoid  effects 
are  ignored. 
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8.  CONCLUSIONS 


With  the  exception  of  the  issues  noted  below,  the  second  phase  of  this  research  is 
complete.  The  formulation  described  in  Kwon  [1]  has  been  implemented,  with  some 
modifications,  into  a  dynamic  finite  element  analysis  program  and  verified.  Damage 
constitutive  equations  have  been  incorporated,  as  well  as  Gurson’s  void  model.  The  effect 
of  void  growth  and  nucleation  on  the  model  has  been  demonstrated,  but  further  verification 
of  the  results  is  still  required. 

The  values  for  the  constants  qb  q2,  and  q3  used  in  Gurson’s  model  are  determined 
by  the  strain  hardening  exponent,  as  described  in  Tvergaard  [4],  A  method  for  determining 
appropriate  values  for  a  given  piecewise  linear  model  needs  to  be  defined.  Incorporating  such 
a  method  into  the  implementation  will  help  simplify  problem  definition,  and  reduce  the  number 
of  constants  the  user  must  provide. 

The  constants  for  void  nucleation  (fN,  eN,  and  sN),  however,  are  dependent  on 
material  type  and  manufacturing  process.  Typical  values  for  a  variety  of  materials  should  be 
provided,  along  with  recommended  default  values. 

The  robustness  of  the  current  implementation  of  the  damage  constitutive  equations 
has  not  been  verified.  Although  this  shell  element  has  been  used  to  analyze  over  20  different 
problems,  further  testing  over  a  wider  range  of  input  conditions  is  needed  to  ensure 
consistency.  Comparing  the  analysis  results  to  experimental  data,  rather  than  the  results  of 
other  models,  should  also  be  accomplished  to  verify  accuracy  and  reliability. 

Work  can  begin  on  the  subsequent  phase  of  this  research,  incorporating  the  element 
into  the  DYSMAS  and/or  DYNA3D  program,  concurrently  with  work  on  the  refinements 
mentioned  above.  The  current  implementation  is  matched  to  a  general  purpose  finite  element 
analysis  program.  All  applicable  portions  of  the  code  will  be  extracted  and  modified  as 
necessary  to  create  a  compatible  module  for  the  target  program. 
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